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Abstract 

In this paper we study a parametric class of stochastic processes 
to model both fast and slow anomalous diffusion. This class, called 
generalized grey Brownian motion (ggBm), is made up of self-similar 
with stationary increments processes (_ff-sssi) and depends on two real 
parameters a £ (0,2) and /3 £ (0,1]. It includes fractional Brownian 
motion when a G (0, 2) and /? = 1, and time-fractional diffusion stochastic 
processes when a = /3 £ (0, 1). The latters have marginal probability 
density function governed by time-fractional diffusion equations of order 
/3. The ggBm is defined through the explicit construction of the underlying 
probability space. However, in this paper we show that it is possible to 
define it in an unspecified probability space. For this purpose, we write 
down explicitly all the finite dimensional probability density functions. 
Moreover, we provide different ggBm characterizations. The role of the 
M- Wright function, which is related to the fundamental solution of the 
time-fractional diffusion equation, emerges as a natural generalization 
of the Gaussian distribution. Furthermore, we show that ggBm can 
be represented in terms of the product of a random variable, which is 
related to the M- Wright function, and an independent fractional Brownian 
motion. This representation highlights the ff-sssi nature of the ggBm and 
provides a way to study and simulate the trajectories. For this purpose, we 
developed a random walk model based on a finite difference approximation 
of a partial integro-differenital equation of fractional type. 

PACS numbers: 02.50Ey, 05.40.Fb, 02.30.Rz, 02.30.Gp 
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1 Introduction 



Diffusive processes are generally classified as normal or anomalous if their 
variance grows linearly in time or not, respectively. Furthermore, the normal 
diffusion is associated to the Gaussian probability density function (PDF) for 
particle positions. 

Several physical phenomena show anomalous diffusion. They range from 
dispersion in complex plasmas |1] to self-diffusion of surfactant molecules j^], 
or from light in a cold atomic cloud [S] to donor-acceptor electron pair within 
a protein [J, to mention only some of the more recent experimental evidences. 
Such anomalous behaviours cover the full range of anomalous diffusion, e.g. 
from slow diffusion ii,, when the variance grows slower than linear, to fast 
diffusion [1] [51 [31 [7j , when the variance grows faster than linear. In order to 
model with a unique mathematical framework both slow and fast anomalous 
diffusion, a class of stochastic processes is here introduced and analyzed. We 
would like this work to be a first step towards a comprehensive description 
of all dispersive mechanisms. Moreover, the PDF of particle positions of this 
class turns out to be related to the M— Wright function, which is a natural 
generalization of the Gaussian density. 

A Gaussian anomalous diffusion can be obtained from a standard diffusion 
equation with time dependent diffusivity. The latter is mainly based on the 
empirical flux-gradient relation and, for this reason, it is considered a simple 
particular case. 

Anomalous diffusion processes can also be obtained as Gaussian processes 
with time subordination (see Remark 1). As a consequence, the particle 
density is not Gaussian. This fact is often seen as the origins and the physical 
interpretation of anomalous diffusion. In fact, let f{x, t) be the density function 
of a diffusive process and G(a;, t) a standard Gaussian density function. Namely, 
with t > 0, 

G(a;,i) = ^^exp f-^V a^t) ■= [ x^G{x,t) dx ^ 2t . (1) 

Let <Pi3{T,t) = t~'^(l){Tt-'^), with T > 0, i > and /3 > 0, be the 
marginal probability density function of a self-similar stochastic process, which is 
interpreted as a randomized operational time. Therefore, in agreement with the 
monotonic growing of time, such a process is required to be a non-negative non- 
decreasing random process. Hereinafter, it is called marginal probability density 
function the one-dimensional PDF of a certain stochastic process (e.g. the one- 
point and one-time PDF of particle positionj^. Furthermore, we remember that 
a process X(t), t > 0, is said to be self-similar with self-similarity exponent H 
if, for all a > 0, the processes X{at), t > 0, and a^X{t), t > 0, have the same 
finite-dimensional distributions. We also suppose that (/) has moments of any 
order. Then, if the subordination formula 

r+oc 

f{x,t) = / Gix,T)ipf3iT,t)dT (2) 

Jo 

^This notation is a physical analogue of probability theory notation. In fact, this marginal 
density corresponds to the result of the integration over an n-dimensional joint density (e.g. 
the n-points and n-times density associated to an n-steps particle trajectory) 



2 



holds, f{x, t) is the marginal density function of an anomalous diffusion process. 
In fact, 



f+oo r+oo ( /'+00 

oo 



hCX2 ( + 



a][t) = / x^f{x,t)dx^ I x'^ \ I G{x,T)ipi3{T,t)dT\ dx 

x^G(x,T)dx^ Lpp{T,t)dT 

= 2 / T^p{T,t)dT = Dtf" , 

Jo 



where we set 

/•+00 

D = 2 cmdc. 

Jo 

which is finite by hypothesis. 

We observe that is desirable having a random time process that ior P — 1 
gives a Gaussian process with a linear growing variance in time. Thus, from [21 
we require that </?i(t, t) ^ 5{t — t). 

A ready-made example is given by the Af-function (see the appendix) , which 
is related to the fundamental solution of the so called time- fractional diffusion 
equation of order (3 (see [U |9l [TOl [TTl [12] ) which, in integral form, reads: 

u{x,t) = uoix) + -— {t-sf ^dxxu{x,s)ds , t>0, (3) 

r(/3) Jo 

with uo{x) — u{x,0). 

In fact, suppose that the marginal density function of the random time 
process is 

(^(r, t) = t~PM0{Tt-P) = Mp{T, t) , 
with < /? < 1. With this choice, [D becomes 

p+ca 

fix,t) = / G{x,T)Ml3iT,t)dT, 

Jo 

which, by using ((15)) and (jlTj) . becomes, 

fix, t) = l M,f2{\xlT)M0{T, t)dT=^ M0/2{\xlt) . (4) 

That is just the fundamental solution of Q. Moreover, when [3 = 1, [3] becomes 
a standard diffusion equation and indeed, using ()44p . one finds 

p + OO 

f{x,t)= G{x,T)Ml{T,t)dT ^ G{x,t) . 



Example 1. Consider a Brownian motion B(t), t >Q, such that E{B{1)'^) = 2 
(we call it standard Brownian motion), and consider a random time l/3{t), t > 0, 
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defined by the local tim(|l in zero at time i of a d = 2(1 — /3)-diniensional 
Bessel process, with < d < 2 (see [131 [H] and [13 [TB]). Furthermore, let Ijs 
be independent oi B. It is know that Z/j is a self-similar process with scaling 
parameter H — (3. Therefore, the subordinated process 

Dp{t)^B{lp{t)), t>0, (5) 

is a model for slow anomalous diffusion and its marginal probability density 
function is the fundamental solutions of the time fractional diffusion equation 
of order < /3 < 1 (see also [T7])- Actually, the local time lj3{t) is defined up 
to a multiplicative constant (see [IH|)- Here, we suppose that lp{t), t > 0, is 
defined such that its marginal density function is M.ii{x, t), x,t > 0. 

Example 2. Consider again a standard Brownian motion B{t), t > 0. Another 
possible choice of an independent random time process lf3(t), for which the 
subordinated process B{lfj{t)) has still marginal density governed by the time- 
fractional diffusion equation of order /3, is the inverse of the totally skewed 
strictly /3-stable process, as founded in the context of Continuous Time Random 
Walk (CTRW) by Meerschaert et al. HH] (see also [101 HI HI 121 IM] ) ■ 

Remark 1. The stochastic interpretation through subordinated processes turns 
out to be very natural. A subordinated process is defined as Y{t) = X{l{t)), t > 
0, where X{t) is a Markovian diffusion and l(t) is a (non-negative) random time 
process independent of X{t) (see [25l|26]). So that, Y{t) has a direct physical 
interpretation. For instance, X{t) can be interpreted as the state of a system 
at time t, while l{t) can be interpreted as the "effective activity" up to time t. 
In this way, even if the process X{t) is Markovian, the resulting subordinated 
process Y{t) is in general non-Markovian and the non-local memory effects are 
attributable to the random time process l{t) and to its evolution, which is in 
general non-local in time. 

Examples 1 and 2 provide two different stochastic processes with the same 
marginal density function Q. Indeed, it is important to remark that, starting 
from a master equation which describes the time evolution of a probability 
density function f{x,t), it is always possible to define an equivalence class of 
stochastic processes with the same marginal density function f{x,t). The two 
examples above represent only particular cases of subordinated type processes. 
However, also processes which are not of subordinated type can serve as models 
for anomalous diffusion described by time-fractional diffusion equations (see 
[27]). It is clear that additional requirements may be stated in order to fix the 
probabilistic model. 

Because < /? < 1, [31 depicts a system with a slow-anomalous diffusion 
behavior. In order to study the full range (slow and fast) of anomalous diffusion, 
we introduce a suitable time-stretching g{t) = t"/^ , < /3 < 1 and a > 0. Let 
f{x, t) be a solution of ([3]), then the function faj3{x, t) — f{x, t"'^^) is a solution 

^Heuristically, the local time l{t,x) of a diffusion process characterizes the "time spent 
by the process at a given level x up to time t" (we shortly write l{t) if x = 0). For 
instance, in the case of Brownian motion, the local time can be written as l{t,x) = 
1 /■* 

lim — / lix-e x+e]i^(^))'t'^i where Ir^ is the indicator function of the inteval. 
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of the stretched time-fractional diffusion equation 

u{x,t) = uo{x) + :^- sf^ ^ [tf^ -s^j —u{x,s)ds, (6) 

with the same initial condition. Then, the fundamental solution of[H]is u{x, t) — 
\Mfi/2{\x\,t°'/P) and defines a self-similar PDF of order H = a/2. That is, 

u(a;,<) = ^Af/3/2(|x|i-"/'), xeR. (7) 

The diffusion is slow when a < 1, standard when a = \ and fast when a > 1. 
We observe that when [3—1, u{x,t) is a "stretched" Gaussian density 

u{x,t)^*-^My,{\x\/t'^/') = ^L=exp(^-^Y t>0. (8) 

Moreover, in the case a — (3, < /3 < 1, the non-Gaussian probability density 
u{x,t) = ^Mij/2{\x\,t) is recovered, i.e. the fundamental solutions of [31 The 
diffusion is always slow and becomes standard when /3 — > 1. Finally, in the 
general case < /9 < 1 and a > 0, we have a non-Gaussian full-ranged 
anomalous diffusion. 



Our main goal is to develop stochastic processes that serve as models for 
the anomalous diffusion described by this class of equations To do this we 
require some constraints. Let X{t), t > 0, be a self-similar stochastic process 
with scaling parameter H = a/2 and marginal probability density function 
defined by [71 We have already observed that there is a whole equivalence class 
of such a stochastic processes. For instance, looking at Examples 1, 2, one 
could immediately take X{t) = B{lp{t°'/^)), t > 0, with a suitable choice of the 
independent random time Ifsit). In order to choose a specific model, we add 
the requirement that the process X(t) be also a stationary increments process. 
Namely, we require that the process X{t) be 77-sssi (self-similar of order H 
with stationary increments process), with H = a/2. 

Remark 2 The latter requirement forces the a-parameter to be in the range 
< a < 2 [28 . Moreover, it automatically excludes the subordinated 
processes of Examples 1, 2, which have in general not stationary increments. 
Summarizing, we ask that the stochastic process X{t), t > 0, satisfies the 
following requirements: let < /? < 1 and < a < 2, then 

i. X{t) is self-similar with index H — a/2. 

l-a/2 

ii. X{t) has marginal density function Ja.p{x,t) — — ^ — Mpi2{\x\t^°'/'^) (see 
0). 

iii. X{t) has stationary increments. 

In [57] the Authors shown that a stochastic process which satisfies all the 
above properties is the so called generalized grey Brownian motion (ggBm) 
Ba,i3{t), t > 0, [5S1 130]. It represents a generalization of Brownian motion 
(Bm) and fractional Brownian motion (fBm) as well. Moreover, it serves as 



5 



stochastic model for [51 So that, in this paper, we will focus on the study 
of this process. Remark 3. Because of the stationarity of the increments, 

the anomalous-diffusion appears deeply related to the long-range dependence 
characterization of Ba.p{t). We remember that an H-sssi process has long- 
range dependence (or long-memory) if 1/2 < _ff < 1. This means that the 
discrete time process of its increments exhibits long-range correlation. That is, 
the increments autocorrelation function 7(fc) tends to zero with a power law 
as k goes to infinity and in such a manner that it does not result integrable 
PSI [5T1 [22]. Therefore, when < a < 1 the diffusion is slow and the process 
has short-memory. While, when 1 < a < 2 the diffusion is fast and the process 
has long-memory. 

The rest of the paper is organized as follows: in the next section we briefly 
introduce the mathematical definition of generalized grey Brownian motion. In 
the third section we characterize the ggBm through the study of its finite- 
dimensional probability density functions. While, the last two sections are 
devoted to trajectory simulations and final remarks. 



2 The generalized grey Brownian motion 

The generalized grey noise space is the probability space (5'(M),S, ^ia,p), where 
5'(M) is the space of tempered distribution defined on R, B is the Borel's a- 
algebra generated by the weak topology on 5'(R) and jia.p is the so called 
generalized grey noise measure. The measure /iQ_/3 satisfies 

e<^^^U^i^,p{Lo) = Ep{-m\l) , ee5(R), 0</?<l, 0<«<2, 

R) 

(9) 

where (•, •) is the canonical bilinear pairing between 5(R) and 5'(R) (see [27] ) 
and Efj{t) is the Mittag-Lefiler function of order (3 

CO „ 

Moreover, 



l^T{l + a) sin -a / dx\x\^-''\f{x)\'' , / G 5(R) , (11) 



with 

J{x) = ^ f dye^yfiv) , / e 

The range < /3 < 1 ensures the complete monotonicity of the function Ei3{—t), 
t>Q (see W R Schneider [33|), as required bylSJ while the range < a < 2 is 
chosen in order to have ||l[a,6)||Q < oo, where l[a.fc) is the indicator function of 
the interval [a, 5). In fact, in this case one has 

\\\a,b)\\l^{b-ar , 0<a<2, Q<a<b. (12) 

It is possible to show (see [27l [SO] [29]) that for each ^ > the real random 
variable 

X.,Ml[o,*))(-)-(-,l[o,t)) (13) 
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is defined almost everywhere on 5'(R). Moreover, it follows from ^ that it 
belongs to L^{S' (R), fia,p) and 

2 

i?(X„,^(l[o,))^) = j^^^r. 
The generalized grey Brownian motion is then defined as the process 

Ba,/3(t) =^a,,3(l[0,t)), * > 0. (14) 

The Bajj{t) marginal density function, indicated with fa^p{x,t), is the 

fundamental solution of [H Namely, fa,p{x,t) — * r[ M^ji{\x\t'°'l'^) (see 
Remark 5). Moreover, the linearity of definition (I14p can be used to show 
many of the fundamentals properties of Ba,f3(t). For instance, Ba^p{t) turns 
out to be i7-sssi with H = a/2. Furthermore, one can calculate characteristic 
functions. For example in the one dimensional case, for any real y and t > 0, 

By using Equations ^ and (fT^ . one has 

LAv,t) = M-y^\\\o.t)\\l) = M-y^n- 

In the multidimensional case, given a sequence of real numbers {9i, 82, ■ ■ ■ , 6n}, 
for any collection {Ba,[3{ti), . . . , Ba^pitn)} with < ti < t2 < • • • < in, using 
linearity again, one can show that P7j 

E I expUy^e.B^Ati)) l^Epl -r(i + /3)i V e,e,j^.Au,t,) I , (15) 



where 



7o,/^(i,^)^ ^^^^^^ ft« + ,"-|t-g|") , t,s>0, (16) 
is the autocovariance matrix of Ba,fj{t). 

Remark 4 From fT5l it follows that, with /? fixed, Ba^p{t) is defined only by 
its covariance structure. In other words, the ggBm, which is not Gaussian in 
general, is an example of a process defined only through its first and second 
moments, which is a property of Gaussian processes indeed. 



3 Characterization of the ggBm 

We want now to characterize the ggBm through its finite dimensional structure. 
From 1151 we know that all the ggBm finite dimensional probability density 
functions are defined only by their autocovariance matrix. The following 
proposition holds 

Proposition 1 

Let Ba^i3 be a ggBm, then for any collection {Bajj{ti), . . . , Ba,p{tn)}, the joint 
probability density function is given by 

{17) 
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with 

1/2 



f = 2r(l y ^^^c.,^3'\U,t 



Proof; in order to show 1171 we calculate its n-dimensional Fourier transform 
and we find that it is equal to pS]) . We have 

/ exp i OjXj fa.fsixi, ... ,Xn; la,l3) d" X = 

, ^^""^ / -^Mb(t) f exp I i | M1/2 ( d'^xdr . 



We remember that Mi/2{r) ~ -^e ^ thus we get 

I -^^p^^) exp J2 ^j^jj ^ 

Vr(l + /3)"det7a,^ "^l^ 2rr(l + /3) j 

We make the change of variables x — r(l + (iY^^T^^^y, with x, y e M", and we 
get 

M^r) exp ^*r(l + ^^y^. ^ x 



y/det -fa, (3 



i.3 = l 



Mp[r) exp |^-r(l + P)t |] ^!2:!i(|lM^ ^ rf^ = e-^^AUr)dr = , 

where s = T{1 + (3) J27j=i (^i(^jla,0{ti,tj)/2 and we have used[43l □ 

Applying the Kolmogorov extension theorem, the above proposition allows 
us to define the ggBm in an unspecified probability space. In fact, given a 
probability space (fi, JF, P), the following proposition characterizes the ggBm 

Proposition 2 Let X{t), t > 0, be a stochastic process defined in a certain 
probability space {Cl,J-,P), such that 

i) X{t) has covariance matrix indicated by 7^.^ and finite- dimensional 
distributions defined by \17\ 
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n) E(X^{t)) 



2 



r /or < /? < 1 and{)<a<2, 



r(i+/3) 



Hi) X{t) has stationary increments, 
then X{t), t>0, is a generalized grey Brownian motion. 

In fact condition ii) together with condition iii) imply that ja.f3 niust be the 
ggBm autocovariance matrix (|16p . 

Remark 5 Using (|^T|l , for n — 1, [T7] reduces to 



This means that the ggBm marginal density function is indeed the fundamental 
solution of [6l 

Remark 6 Because ioi P — 1 

Mi(t) = (5(t- 1), 

then, putting "/a,i = 7q, we have that 1171 reduces to the Gaussian distribution 
of the fractional Brownian motion. That is, 

fa,i{xi,X2, ■ ■ . ,a;„;7a,i) = ^2det7 ^^^^^ I I ^ ^ ^«7a ^(^i: 1 j ■ 
We have the following corollary. 

Corollary 1 Let X{t), t > 0, be a stochastic process defined in a certain 
probability space {Q,J-,P). Let H — a/2 with < a < 2 and suppose that 
E{X(1)'^) = 2/r(l + P). The following statements are equivalent 

i) X is H-sssi with finite- dimensional distribution defined by {11^ , 

ii) X is a generalized grey Brownian motion with scaling exponent a/2 and 
"fractional order" parameter (3, 

iii) X has zero mean, covariance function Ja,i3{t,s), t,s>0, defined by fT^) 
and finite dimensional distribution defined by ( |j7| ). 

3.1 Representation of ggBm 

Up to now, we have seen that the ggBm Ba^p{t), t > 0, is an _ff-sssi process, 
which generalizes Gaussian processes (it is indeed Gaussian when [3 — \) and 
is defined only by its autocovariance structure. These properties make us think 
that Ba^i3{t) may be equivalent to a process Ki}Xa{t), t > 0, where Xa{t) is 
a Gaussian process and is a suitable chosen independent random variable. 
Indeed, the following proposition holds 



fa,l3{x,t) 




(19) 
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Proposition 3 Let Ba^pit), t>Q, he a ggBm, then 

Sa,/3W = TAa^aW , i>0, 0</3<l, 0<a<2, (20) 

where = denotes the equality of the finite dimensional distribution, Xa{t) is 
a standard fBm and is an independent non-negative random variable with 
probability density function Mp{T), r > 0. 

In fact, after some manipulation, [TH] can be written as follows 

/ / exp f i V djyxj ] 2yM^(y2)ijL4 exp f - V x,-f-\t^,tj)xj /2 ] dyd^'x 



E |^exp(i^0,V^^afe)) 

Example 3 A possible choice of is the random variable lp{l), where lp{t), 

i > 0, is the random time process of Example 1 or Example 2. 

Remark 7 Proposition [3] highlights the iJ-sssi nature of the ggBm. Moreover, 
for P — 1 from El follows that Li = 1 a.s., thus we recover the fractional 
Brownian motion of order H — a/2. 

Representation (j20p is very interesting. In fact, a number of question, in 
particularly those related to the distribution properties of Ba^f3{t), can be 
reduced to questions concerning the fBm (t) , which are easier since Xa (t) is a 
Gaussian process. For instance, the Holder continuity of the Ba^f3{t) trajectories 
follows immediately from those of Xa (t) 



EQX^it) - X^is)\P) = cp\t - s 



pa/2 



Moreover, this factorization is indeed suitable for path-simulation (see next 
section). 

Remark 8 From (|20p . it is clear that the Brownian motion {Bi^i{t), t > 0) is 
the only one process of the ggBm class with independent increments. 



4 Path simulation 

In the previous section we have shown that the ggBm could be represented by 
the process 

Ba,j3{t) = y/L^Xa{t), t>0, 0<(3<1, 0<a<2, 

where is a suitable chosen random variable independent of the standard 
fBm Xa{t). Clearly, to simulate ggBm trajectories we first need a method to 
generate the random variable Lp. 
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4.1 The time fractional drift equation 



In order to generate the random variable Lfj with probabihty density function 
Mjs^r), we consider the so called time-fractional forward drift equation, which 
in integral form reads 

1 /■* d 

u(x,t)^uo(x) — (t-s)'^-^—u{x,s)ds, xeR. t>0, < (3 < 1 . 

r(/3) 7o ox 

(21) 

The fundamental solution of[^is ^TE[ [TC] 

u{x,t) = Mf3{x,t) , x,t>0. (22) 

This function can be interpreted as the marginal density function of a non- 
negative self-similar stochastic process with scaling parameter H = [3 (see 
Examples 1, 2). 

Remark 9 The name "drift equation" refers to the fact that when (3 — ll^turns 
out to be the one dimensional (forward) drift equation dtu{x,t) = —dxu{x,t), 
whose fundamental solution is d{x ^ t). 

Remark 10 When /3 = 1 we recover Mi{x, t) — S{x — t) (see HH)). 

We write (^11 in terms of fractional derivative of order 13. Let us introduce 
the Caputo-Dzherbashyan derivative 

where J" is the Riemann-Liouville fractional integral of order < a < 1 such 
that, for a = 0, is the identity operator. Then, the corresponding Cauchy 
problem of[5T]can be written 

*D'^u{x,t) = ~da:U{x,t) , ^24) 

u{x,0) = uo{x) ~ 6{x) , 

with a; e R, t > and < /? < 1. 

Using a random walk model, one can simulate a discrete time random process 
Lp{t), t > 0, governed by the time- fractional forward drift equation (see 
[341 135] ). In this way, for each run, the random variable Lfj{l) has the required 
distribution u(x, 1) = M0{x). The random walk construction follows two steps: 

• the Griinwald-Letnikov discretization of Caputo-Dzherbashyan derivative, 

• the interpretation of the corresponding finite difference scheme as a 
random walk scheme. 



4.2 Finite difference schemes 

In order to define the finite difference model, we write the Cauchy problem 
in a finite domain 

*L)fM(x,t) = -^w(a;,t) , (x, t) e O = [-a, a] x [0, 1] , a>0, 
u{x,Q) = Ui^{x) = 5{x) , (25) 
u(-a,t) = $i(t) , u(a,t) = $2(i) , t>Q. 
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Let iV, M be positive integers. Then, we introduce a bi-dimensional lattice 
Ss^st' ^ {USx,nSt), {j,n) e Z2M x Zjv} , 

contained on Q, with 6x = 2a/{2M - 1) and 6t = 1/{N - 1). The lattice 
elements are indicated with 

{xj,tn) ^ {jSx,nSt) , j = 0, 1, . . . ,2Af - 1 , n = 0, 1, . . . iV - 1 . 

Let u : ^ M be a function defined on H.. We indicate with it" = u{xj,tn) the 
restriction of u to Gs^fst' evaluated in {xj,tn). 

The time-fractional forward drift equation is then replaced by the finite 
difference equation 



where 



^^f-"--^^^^> (26) 



fc=0 



(27) 



/3\ r(/3 + 1) 



^fcy r(/c + i)r(/? - /c + 1) ' 

is the so called forward Griinwald-Letnikov scheme for the Caputo- 
Dzherbashyan derivative (|23p . Using the "empty sum" convention 

E • = , if g < p , 

k—p 

for any n > 0, we obtain the explicit equation 

= E(-i)'^ (f ) + tyr' {i) -r-' + A^K-i - (28) 

where ^ = St^ /6x. Equation (|28p can be written in the following noteworthy 
form 



uj+i = 6„u^" + E Cfcii;'+^-'^ + - uj), (29) 

where we have defined 



(30) 



More precisely, the explicit scheme reads 

U0=$l(tn), W2A/~1 = ^2(tn), n > 0, 

u] = (1 - /i + fiL)u^j, 0<j< 2M - 1 

[ci- fl + II^U"^ + C2vJf^ + ...+ CnU] + bnU"j, Ti > 0, < J < 2Af - 1, 
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where L is the "lowering" operator Lfj = fj-i- 

Remark 11 When /3 = 1 all the coefficients Cfc and 6„ vanish except 60 = ci = 1. 
So that, we recover the finite difference approximation of the (forward) drift 
equation. 

It order to write an impUcit scheme, we have to use backward approximations 
for the time-fractional derivative. Therefore, for any n > 0, we obtain 

n 

U]+' + I,{U';+' - U]+1) = hnvPj + J2 CkU]+'-'. (31) 

k=l 

Namely, 

?i[J = $i(i„), u^j^_^^ ^2{t„), n>0 
(l + /u-/xX)u] =1*0, 0<i<2M-l, 

(1 + M - IJ^L)u]+^ = ciu'^ + C2uJ-^ + . . . + Cnu] + n > 0, < j < 2M - 1. 
The above equation can be rewritten in matrix notation 

Aiju]+^ = CiUf + C2<-' + . . . + CnUj + bnVpi + V-f+S n>l. ^^^^ 

A is the following 2M x 2M matrix, divided in four M x M blocks 

A = 



Ai 





A2 


A ^ 



A = 



1 



















... 




1 + n 
















... 





-M 


1 + fj. 














... 










-jJ- 1 + A* 
























1 + /i 








... 














-M 


1 + At 





... 



















1 + M 


... 






















1 



(33) 

Moreover, for any n > 0, V" is a suitable vector which takes into account for 
the boundary terms. 

Remark 12 Because A is lower diagonal, A~^ is 





~ 


-A-iAaAf' 





(34) 



As usual, the explicit scheme is subjected to a stability condition, while te 
implicit scheme is always stable. For example, if we take n < (3, namely 

5x > 6t^l(3, (35) 
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then the exphcit scheme results indeed stable and preserves non-negativity as 
well. This means that if we suppose u° > for any < j < 2Af — 1, then 
-u" > for any n > and < j < 2M — 1. Actually, this is crucial because 
we interpret the {u"} as sojourn probabilities. In order to show this, it is 
convenient to write equations (j29|) and (I3ip in the Fourier domain. Namely, 
we apply the discrete Fourier transform with respect to the "space" index j to 
both sides of and (pij) . We remember that, given a collection of complex 
numbers {xj, j = 0, l,---,2Af — then its discrete Fourier transform is 
usually defined as: 

J'd{xi)k:=Xk= 2:,e-»2.jV2M^ fc = 0, • • • , 2M - 1. 

3=0 

One can show that, for any real number a, one has: 

Thus, heuristically, the effect of applying the Fourier transform to our finite 
difference equations is just to "line up" the points in the Fourier space. In fact, 
for any n > 0, we get: 

u^~^^ — ie{k)u^ + (• • explicit case, 

ii{k) — Pu^ + (• ■ ■)k, imphcit case, 

where Ce(fc) = — + A"3~*'^'^^*^) and^i(fc) — (1 + ^ - /xe~*'^'^/^^)~^ correspond 
to the so called amplification factors. Then, heuristically, the schemes are stable 
if I < 1 for any k. Namely, one has to require: 

max ((/? — + + 2^(/3 — fi) cos{nk/M)) < 1, explicit case, 
min ((1 + yu)^ + /i^ — 2/_t(l + ^) cos(7rfc/Af )) > 1, imphcit case. 

k 

Clearly, while the second one is always satisfied, the first condition is indeed 
true if dSni) holdfl 

4.3 Random walk models 

We observe that, for < /? < 1, 

oo 

^ Cfc = 1, 1 > /? = ci > C2 > . . . ^ . (38) 

k=l 

Moreover, 

oo m oo 



fc=l fc=l k=m+l 

1 = 6o > 6i > 62 > • • • -> . 



(39) 



Thus, the coefficients Ck and 6„ are a sequence of positive numbers, which do 
not exceed unity and decrease strictly monotonically to zero. 



^Actually, we are considering {xj}, for any j G Z-f, as a periodic sequence, such that 
Xj+2M = Xj for any non- negative integer j. 

^This follows from the fact that if (/3 — /x) is non-negative, the maximum is reached when 
the cosine equals +1, then one has < 1, which is in fact true by hypothesis. 
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4.3.1 Explicit random walk 



In order to build a random walk model, we consider for first the explicit scheme. 
Omitting the bomidary terms, we have 

r u] = (i-/i)MO + /iw°_i, 

\ U]+^ = (Ci - + + C2U]-^ + ...+ Cnu] + bnU'l , n>l. 

We consider a walker which starts in zero at time zero, namely x(t = 0) = 0. 
We interpret the u" as the probability of sojourn in Xj = jSx at time t„ — nSt. 
Then, we indicate with x{tn) the position of the particle at time i„. 

At time ti = St the walker could be at the position a;(l) = xi with probability 
fi (that is the probability to come from one space-step behind) or in — 
xq = with probability 1 — /i (that is the probability to remain in the starting 
position). 

From Equations psp and (|39p. it is clear that the parameters ci, C2, . . . c„, 6„ 
can be interpreted as probabilities. Then, the position at time tn+i is 
determined as follows. We define a partition of events {Eci , , ■ ■ ■ , Ec„ , i?6„ }, 
with P{Ec^) — Ck, P{Eb^) = bn, n > 1 and such that: 

• i?ci —^the particle starts in the previous position x(tn) and jumps in 
x{tn) + 5x with probability ji or stays in x{tn) with probability 1 — 

• i?cfc =|t/ie particle backs to the position x(t„+i_fe)|. 

• Eb^ —^the particle backs to the initial position x(to)|- 

4.3.2 Implicit random vi^alk 

Consider the implicit case 
r = A-lu^ 

\ = [ciu" + C2w""^ + . . . + c„ui + b„u°] , 

where A is given by (|34p . The probability interpretation of the parameters 
and bn is still valid. In this case, however, we must use the transpose matrix P 
of to define the transition probabilities. Indeed, the M x M matrix 
propagates in the positive semi-axis. Moreover, from ()33|) . it can be shown that 
P defines a transition matrix (all the elements are positive numbers less than 
one and all the rows sum to one, see example below). 

Example 4 In the four dimensional case the matrix P is 

/i/(i + M) m/(i + m)' ^^V{l + ^^? + 
1/(1 + A^) m/(i + m)' mV(i + m)' 

1/(1 + A^) A^/(l + M) 

\ 1 / 

In the implicit case the random walk is defined as follows. Let the particle 
start in zero, at the first step it could jump up to M — 1 steps ahead with 
probabilities defined by the first P row, then we have the following partition of 
events: 
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• =|i/ie particle starts in the previous position x{tn) for instance Xj. 
Then, it could jump up to M — j — \ steps ahead with probabilities defined 
by the j + 1-th P row^\ 

• i?Cfc —^the particle backs to the position a;(t„+i_fc) for instance Xj. Then 
it could jump up to M ~ j I steps ahead with probabilities defined by the 
j + l-th P row^; 

• Eh,^ ^^the particle backs to the initial position x{tQ) then could jump up 
to M — 1 steps ahead with probabilities defined by the first P row^ . 

Remark 13 The implicit method is slower than the explicit one. However, 
we observe that, because of the stability constraint ((35|) . the implicit scheme is 
advisable for small (3. Indeed, we have 5x ~ St^ /P, with < (3 < 1, and we are 
forced to raise the time steps in order to improve "spatial" resolution. 

4.4 ggBm trajectories 

Along this section, we have provided a method to generate the random variable 
Lp and the simulation results are shown in Figures 1-3. In particular 

• Figure 1 shows a random-walk simulation with (3 — 0.4. In this case 
we used an implicit random- walk scheme. Moreover, we compared the 
histogram evaluated over N = 10000 simulations and the density function 

M2/5ix), x>0. 

• Figure 2 shows a random walk when /3 — 0.5. Because Mi/2ix) = 
_^g-a; 2: > 0, in this case L1/2 — \Z\ where Z is a Gaussian random 
variable. 

• Figure 3 shows a random walk with P = 0.8. We used an explicit random- 
walk scheme. Then, we compared the histogram with the corresponding 
probability density function M4/5(a;), a; > 0. 

In all the studied cases, we have found a good agreement between the 
histograms and the theoretical density functions. 

At this point, in order to obtain examples of the Ba^p{t) — ^ LpXa{t) 
trajectories, we just have to simulate the fractional Brownian motion A"„(t). 
For this purpose we have used an exact Cholesky method (see [36j ) . 

Path simulations of Bp^ij{t) (shortly Bp) and B2-i3^p{t), with (3 = 1/2 are 
shown in Figures 4-9. The first process provides an example of stochastic model 
for slow-diffusion (short- memory), the second provides a stochastic model for 
fast-diffusion (long- memory) . 

• Figure 4 shows some typical paths. In the bottom panel, we present the 
corresponding increment process. Namely, 

Zp{tk)^Bp{tk)-Bp{tk-i), tk = kSt, fc= l,2,...,Af- 1. 
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Figure 1: In the top panel, the histogram of L^, which is calculated from 
a sample of = 10000 outcomes, is obtained with an implicit random-walk 
scheme and it is compared with the exact PDF Mis{x), x > 0, with /3 = 0.4. In 
the bottom panels, the random variable Lp (left) and two trajectory examples 
(right) are shown. 
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Figure 2: In the top panel, the histogram for the case f3 — 0.5 is calculated from 
a sample of = 15000 outcomes, which are obtained simulating independent 
Gaussian random variables. The corresponding Lp is shown in the bottom. 



• Figure 5 shows the agreement between simulations and the theoretical 
densities at times t — 1 and t = 2. 
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Figure 3: In the top panel, the histogram of Lp, which is calculated from a 
sample of = 15000 outcomes, it is obtained with an explicit random-walk 
scheme and it is compared with the exact PDF M^(x), x > 0, with /3 = 0.8. In 
the bottom panels, the random variable Lp (left) and many trajectory examples 
(right) are shown. 




t 



Figure 4: Bp{t) trajectories in the case /3 = 0.5 (top panel) for < t < 2. 
The time series of the corresponding stationary noise Zf}{t) is presented in the 
bottom panel. 
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Figure 5: Histograms of = 15000 simulations with (3 = 0.5 and exact marginal 
density at t = 1 and t = 2. 



• Figure 6 presents the plot of the sample variance in logarithmic scale. 
Moreover, we evaluated a linear fitting, which shows a good agreement 
with the theoretical result. 

• Figure 7 shows some typical paths for the long- memory process -B2-/3./3(i). 

• Figure 8 collects the histograms in the case a = 2 — /3 at time t = \ and 
time t — 2. The super-diffusive behavior (that is the rapid increasing of 
the variance in time) is highlighted. 

• Figure 9 presents the plot of the sample variance in logarithmic scale. 
Even in this case, we evaluated a linear fitting. 



5 Concluding remarks 

The marginal probability density function (I19p of the generalized grey Brownian 
motion Ba^f3{t), t > 0, evolves in time according to a "stretched" time- 
fractional diffusion equation of order /? (see[S]). Therefore, the ggBm serves 
as stochastic model for the anomalous diffusion described by these class of 
fractional equations. 

The ggBm is defined canonically (see [T4| in the so called grey noise space 
{S'{R),B, fj.a.,13), where the grey noise measure satisfies However, the 

ggBm is an ff-sssi process of order H = a and Proposition [5] provides a 
characterization of Ba.p{t) notwithstanding the underline probability space. 

There are many other processes which serve as stochastic models for a given 
master equation. In fact, given a master equation for a PDF f{x,t), it is always 
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Figure 6: Sample variance in logarithmic scale and linear fitting {N = 10*). 
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Figure 7: B2-p.i3{t) trajectories in the case (3 = 0.5 (top panel) for < t < 2. 
The corresponding stationary noise time series is presented in the bottom panel. 



possible to define an equivalence class of stochastic processes with the same 
marginal density function /(a;, t). The ggBm defines a subclass {Ba^pit)^ i > 0} 
associated to the non-Markovian equation In this case, the memory effects 
are enclosed in the typical dependence structure of a iJ-sssi process. While, 
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for instance in the case of a subordinated process (Examples 1 and 2), these 
are due to the memory properties of the random time process. The latter are 
preferable because provide a ready-made physical interpretation (Remark 3). 
However, Ba^f3{t) is interesting because of the stationarity of its increments. 
Proposition [3] provides an enlighten representation of ggBm. Thus, the 
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generalized grey Brownian motion turns out to be merely a fractional Brownian 
motion with stochastic variance, that is Ba^p{t) = KpXa{t), t > 0, where 
is a suitable independent random variable (see [5D]). As a final remark, 
we observe that such a process is not ergodic, as, heuristically, follows by 
the multiplication with the random variable A^. This appears also from the 
simulated trajectories. Indeed, it is impossible with a single realization of the 
system Ba,i3{t,u}), t^; G il, to distinguish a ggBm from a fBm with variance 
2A|(a;)t", where A^(ti;) indicates a single realization of the random variable A/3. 

As a further development of the present research, it is interesting to wonder 
if the generalized grey Brownian motion is the only one stationary increment 
process which serves as model for time- fractional diffusion equations like ©. 
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Appendix. The M-function 

Let us define the function Mp{z), < /3 < 1, as follows 

= -E^^r[(/3(fc + l))]sin[^/3(fc + l)] , r>0. (40) 

TT ^ — ' fc! 

fe=0 

The above series defines a transcendental function (entire of order 1/(1 — (3/2)) 
of the Wright type, introduced by Mainardi in [23 US], (see also [HE]). It is 
useful to recall some important properties of the M-function. The best way to 
express these properties is in terms of the function 

which defines a probability density function in t > for any < > and < /3 < 1. 

1. The following convolution formula holds: 

Mu{x,t)= I Mjj{x,T)Mi3{T,t)dT , v = r]l3, x > , (41) 

where v.rj.P ^ (0, 1). 

2. The Laplace transform of Aip{T,t) with respect to t is 

C{Mp{T,t)-t,s) = sl^-^e-^'^ r,s>0. (42) 
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3. The Laplace transform of A^^(t, t) with respect to r is 

£{Mi3{T,t);T,s}^Ep{-st'^), t,s>0, (43) 
where Ep{x) is the Mittag-Leffler function (fTU)) . 

4. The singular limit f3 1 gives 

Mi{T,t) ^ S{t -t) , r,i>0. (44) 

5. Let G{x, t) be the Gaussian density function, then 

G{x,t) - ixi/2(|a;|,t) = ^exp (^-^^ . (45) 
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